rm(list = ls())
options(width=100, device = 'quartz')

###################
## Assemble the Data 
####################
library(foreign)
library(plyr); library(dplyr)
library(coin)
library(rdrobust)
library(rms)
library(rdd)
library(ggplot2)
library(reshape2)
library(xtable)
library(MASS)
library(plm) ## for panel models
library(stargazer) ## for nice latex tables
library(lmtest) ## for coeftest
library(clusterSEs) ## for bootstrapped SEs

myRecode <- function(var, recode.ls, as.factor.result=TRUE, ...) {
  options(useFancyQuotes = FALSE)
  recode.ls <- sapply(recode.ls, function(x) paste(sQuote(x), collapse='='))
  recodes <- paste(recode.ls, collapse=';')
  print(recodes)
  revar <- Recode(var, recodes, as.factor.result=TRUE, ...)
  options(useFancyQuotes = TRUE)
  revar
}
rols <- function (..., data, subset=TRUE, cluster=1:nrow(data),
                  method="huber") {
  mod <- ols(..., data=data[subset, ], x=TRUE)
  robcov(mod, cluster=cluster[subset], method=method)
}
PasteEval <- function(..., print = FALSE) {
  txt <- paste(..., sep = "")
  if (print) { cat("\n", txt, "\n") }
  eval(parse(text = txt), envir = parent.frame())
}
FileName <- function (path=NULL, name=NULL, ext=NULL, replace=FALSE) {
  p <- paste(path, collapse = "")
  n <- paste(name, collapse = "")
  e <- paste(".", ext, sep = "")
  d <- format(Sys.Date(), "%y%m%d") 
  for (l in seq_along(letters)) {
    if (l > 1) old_file_name <- file_name
    file_name <- paste0(p, d, paste(n, letters[l], sep="-"), e)
    if (!any(grepl(file_name, list.files()))) {
      if (replace) file_name <- ifelse(l > 1, old_file_name, file_name)
      break
    }
  }
  cat("\nFile name:", file_name, "\n\n")
  return(file_name)
}

mysw <- function (...) {
  setwd(paste(..., sep = "/"))
}
mypdf <- function(name, ...) {
  date <- format(Sys.Date(), "%y%m%d")
  pdf(paste(paste(name, collapse=""), date, ".pdf", sep=""), ...)
}
tri <- function (x, h, c=0) pmax(0, 1 - abs((x - c) / h))
POtoSouth11 <- function (pos, dta=st.info) {
  dta$South11[match(pos, dta$POAbrv)]
}



mysw(home.dir)

st.info <- read.csv("StateCodes.csv")


load(file= "party_control_data_161121.RData")

data.use <- data %>%
    plyr::arrange(abb, year) %>%
        dplyr::filter(!is.na(abb) & year >= 1936)

table(data.use$year)

################
#### Setup RDD Variables
################

data.use <- data.use %>%
    group_by(
        abb
    ) %>% mutate(
        South11 = POtoSouth11(abb),
        DemGov = gov_party,
        DemPropGov0 = demprop2 - 1/2,
        DemMarginGov = DemPropGov0 * 2,
        DemWinGov = as.integer(DemPropGov0 > 0),
        DemSeatShareHouse = hs_dem_per_2pty / 100,
        DemSeatShareSenate = sen_dem_per_2pty / 100,
        DemControlHouse = hs_dem_control,
        DemControlSenate = sen_dem_control,
        GovLib = Policy,
        GovLibL1 = lag(GovLib, 1),
        GovLibL2 = lag(GovLib, 2),
        GovLibL3 = lag(GovLib, 3),
        GovLibL4 = lag(GovLib, 4),        
        GovLibL12 = (GovLibL1 + GovLibL2) / 2,
        GovLibP1 = lead(GovLib, 1),
        GovLibP2 = lead(GovLib, 2),
        GovLibP3 = lead(GovLib, 3),
        GovLibP4 = lead(GovLib, 4),
        GovLibP12 = (GovLibP1 + GovLibP2) / 2,
        GovLibD1 = GovLibP1 - GovLib, ## delta 1
        GovLibD2 = GovLibP2 - GovLib,  ## delta 2
        GovLibD3 = GovLibP3 - GovLib,  ## delta 3
        GovLibD4 = GovLibP4 - GovLib,  ## delta 4
        GovLibD12 = (GovLibD1 + GovLibD2)/2, ## ave of deltas 1-2
        GovLibDL1 = GovLib - GovLibL1 ## change, t relative to t-1
    ) %>% ungroup()


res.fe <- resid(lm(GovLib ~ abb + factor(year), data=data.use))
res.l1 <- resid(lm(GovLib ~ Policy_L1, data=data.use))
res.l2 <- resid(lm(GovLib ~ Policy_L2, data=data.use))
res.fe.l1 <- resid(lm(GovLib ~ abb + factor(year) + GovLibL1, data=data.use))
res.fe.l2 <- resid(lm(GovLib ~ abb + factor(year) + GovLibL2, data=data.use))

data.use$GovLib.fe0 <- NA
data.use$GovLib.fe0[as.integer(names(res.fe))] <- res.fe
data.use$GovLib.fe1 <- NA
data.use$GovLib.fe1[as.integer(names(res.fe.l1))] <- res.fe.l1
data.use$GovLib.fe2 <- NA
data.use$GovLib.fe2[as.integer(names(res.fe.l2))] <- res.fe.l2
data.use$GovLib.1 <- NA
data.use$GovLib.1[as.integer(names(res.l1))] <- res.l1
data.use$GovLib.2 <- NA
data.use$GovLib.2[as.integer(names(res.l2))] <- res.l2

data.use <- group_by(data.use, abb) %>%
    mutate(GovLib.1L1 = lag(GovLib.1, 1), ## lag of first-difference
           GovLib.1L2 = lag(GovLib.1, 2),
           GovLib.2L1 = lag(GovLib.2, 1), 
           GovLib.2L2 = lag(GovLib.2, 2),
           GovLib.fe0L1 = lag(GovLib.fe0, 1),
           GovLib.fe0L2 = lag(GovLib.fe0, 2),
           GovLib.fe1L1 = lag(GovLib.fe1, 1),
           GovLib.fe1L2 = lag(GovLib.fe1, 2),
           GovLib.fe2L1 = lag(GovLib.fe2, 1),
           GovLib.fe2L2 = lag(GovLib.fe2, 2),
           GovLib.1P1 = lead(GovLib.1, 1), ## lead of first-difference
           GovLib.1P2 = lead(GovLib.1, 2),
           GovLib.2P1 = lead(GovLib.2, 1), 
           GovLib.2P2 = lead(GovLib.2, 2),
           GovLib.fe0P1 = lead(GovLib.fe0, 1),
           GovLib.fe0P2 = lead(GovLib.fe0, 2),
           GovLib.fe1P1 = lead(GovLib.fe1, 1),
           GovLib.fe1P2 = lead(GovLib.fe1, 2),
           GovLib.fe2P1 = lead(GovLib.fe2, 1),
           GovLib.fe2P2 = lead(GovLib.fe2, 2)
           )

tail(dplyr::select(data.use, c(year, GovLib, GovLib.fe0, GovLib.fe0L1, GovLib.fe0L2,
                              GovLib.fe0P1, GovLib.fe0P2)))


#### PARTISAN DIVERGENCE
summary(data.use$StateHouseDivergence_L1)
summary(data.use$StateSenateDivergence_L1)

data.use <- data.use %>% mutate(
    StateLegDivergence =
        (StateHouseDivergence + StateSenateDivergence) / 2,
    StateLegDivergence_L1 =
        (StateHouseDivergence_L1 + StateSenateDivergence_L1) / 2
)

mypdf("Figure7_MassEliteDivergence", width=6, height=6)
(ggplot(data.use %>% filter(
    !is.na(StateLegDivergence)
) %>% group_by(
    abb
) %>% summarise(
    MassPartisanDivergence = mean(MassPartisanDivergence, na.rm=TRUE),
    StateLegDivergence = mean(StateLegDivergence, na.rm=TRUE)
))
 + aes(x=MassPartisanDivergence, y=StateLegDivergence, label=abb)
 + geom_text(alpha=1/2)
 + geom_smooth(method="lm", formula = y ~ splines::bs(x, 3), se=FALSE)
 ## + geom_smooth(method="lm", se=FALSE)
 + labs(x="Partisan Divergence in the Mass Public",
        y="Partisan Divergence in the State Legislature",
        title="Mass and Elite Partisan Divergence, 1993-2014")
 + theme_bw()
 )
dev.off() 
 
